function [spectra, descriptions] = ...
  generateFlowFrequency( flow, features, Fs )
%generateFlowFrequency - Computes the flow frequency.

  if ~exist('Fs', 'var')
    Fs = 25;
  end
  
  % Initialize some variables
  N = size(flow, 1);
  spectra = zeros(length(features)*2, N);
  descriptions = cell(length(features),1);
  
  % Iterate through each feature, transforming each into
  % a frequency spectrum
  addpath utils
  for f=1:length(features)
    i = (f-1)*2 + 1;
    [Y, desc] = generateFlowFeatures(flow, N, features(f));
    spectra(i,:) = fourierTransform(Y(:,1)', N, Fs, false);
    spectra(i+1,:) = fourierTransform(Y(:,2)', N, Fs, false);
    descriptions{f} = desc;
  end
  
end

